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I. 


Introduction 


An  early  and  continuing  goal  of  lidar  research  has  been  to 
devise  an  inversion  method  whereby  profiles  of  optical  parameters 
such  as  attenuation  and  backscatter  coefficients  in  an  inhomogenous 
atmosphere  can  be  quickly  and  accurately  deduced  from  the  return 
signal  of  a  monostatic,  single  wavelength  lidar  system.  This  is  a 
problem  area  where,  as  expressed  by  Collis  and  Russell  in  an  excellent 
review  article  ”...  the  early  promise  of  lidar  has  not  yet  been 
fulfilled”  [1].  Some  of  the  difficulties  encountered  along  the  way 
have  been  due  to  limitations  in  lidar  performance  and  associated 
data  processing  technology,  while  others  follow  from  theoretical 
requirements  and  constraints  peculiar  to  the  inversion  process. 

This  article  addresses  some  aspects  of  the  latter  category  of 
problems,  and  presents  in  particular  a  simple  inversion  method 
based  on  a  new  form  of  a  well-known  analytical  solution. 

II.  Review  of  the  "Slope”  and  "Solution”  Methods  of  Inversion 

For  a  monostatic,  single  wavelength,  pulsed  lidar,  the  assumed 
basic  governing  form  is  the  single-scattering  lidar  equation: 

P(r)  *  Po  p  A  expl-l^'ofr'Jdr'J  (1) 

In  this  equation  P(r)  is  the  instantaneous  received  power  at  time 
t,  Po  the  transmitted  power  at  time  tQ,  c  the  velocity  of  light,  t 
the  pulse  duration,  A  the  effective  system  receiver  area,  r(*c(t-to)/2) 
is  the  range,  and  £(r)  and  o(r)  are  respectively  the  volume  backscatter 
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and  atCenuation  coefficients  of  the  atmosphere.  A  more  convenient 
signal  variable  is  the  logarithmic,  range-adjusted  power,  defined 
as 


S(r)  =  in{r2P(r) ]  (2) 

In  terms  of  S  =  S(r)  and  So  =  S(rQ),  where  r^  is  a  given  constant 
reference  range,  Eq.  (1)  may  be  expressed  in  a  system- independent 
form: 


r 


r 

o 


where  pQ  =  P(rQ). 


The  differential  equation  corresponding  to  Eq.  (3)  is 


dS  1  dg 
dr  P  dr 


20, 


(4) 


a  solution  which  evidently  requires  knowing  or  assuming  a  relation¬ 
ship  between  P  and  a  whenever  dp/dr  t  0.  On  the  other  hand,  if  the 
atmosphere  is  homogeneous  so  that  dp/dr  =  0,  the  attenuation  coef¬ 
ficient  can  be  expressed  directly  in  tens  of  the  signal  slope: 


_  J.  dS 
hom.  ”  2  dr 


(5) 


This  is  the  basis  of  the  slope  method  of  inversion  [2,3],  in  which 
typically  the  slope  of  the  Least  squares  straight  line  fit  to  the 
curve  S  =  S(r)  is  used  as  the  best  estimate  of  dS/dr  over  any 
interval  where  S  itself  appears  to  be  nearly  a  straight  line. 
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Going  a  step  further,  it  has  often  been  assumed  that  since  the 
atmosphere  is  more  likely  to  be  homogeneous  over  small  rather  than 
large  intervals,  by  applying  the  slope  method  to  a  succession  of 
small  intervals  a  reasonable  first  approximation  to  o  =  a(r)  in  a 
notably  inhomogenous  atmosphere  may  also  be  achieved.  From  Eq.  (4) 
it  is  clear  this  amounts  to  a  conjecture  that  generally  p  * J dp/dr |  <<2a, 
at  least  over  most  of  the  S  curve.  Unfortunately,  assumptions  like 
this  appear  not  to  be  well  justified  for  many  situations  of  interest, 
e.g.,  under  conditions  of  dense  cloud,  fog,  smoke,  and  dust.  Even 
under  the  relatively  stable  conditions  prevailing  in  fogs,  significant 
local  heterogeneities  occur.  For  example,  the  spatial  variation  of 
fog  drop  concentrations  is  often  quite  large,  ranging  up  to  two 
orders  of  magnitude  for  certain  size  categories  [4,5].  Such  micro¬ 
structure  variation  along  the  lidar  beam  path  could  easily  lead  to 
relatively  large  fluctuations  in  dp/ dr,  hence  invalidating  local 
application  of  the  slope  method.  The  same  criticism  applies  to  the 
so  called  "ratio”  or  "slice"  method  of  inversion  [6,7],  which  is 
merely  an  extremely  close  variant  of  the  slope  method  as  applied  to 
successive  range  intervals.  (Additional  discussion  on  the  merits 
of  the  slope  and  ratio  methods  is  available  through  the  recent 
articles  of  Kohl  [8,  9]  and  Brown  [10].  Also,  a  theoretical  example 
which  illustrates  the  inadequacy  of  the  slope  method  for  a  case  of 
high  visibility  is  given  below  in  Section  V.) 

Several  observational  and  theoretical  studies  have  been  published 
which  show  that  under  a  wide  range  of  realistic  circumstances  P  and 
a  can  in  fact  be  related  approximately  according  to  a  power  law  of 
the  form 

P  *  const.  ak, 
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(6) 


where  k  depends  on  the  lidar  wavelength  and  various  properties  of 
the  obscuring  aerosol.  Reported  values  of  the  exponent  are  generally 
on  the  interval  0.67SkS1.0  [11-15].  If  such  a  relationship  is 
assumed,  Eq.  (4)  becomes 


dS  =  k  do 
dr  ”  o  dr 


2  a 


(7) 


Although  the  above  ordinary  differential  equation  is  nonlinear, 
it  nevertheless  has  an  elementary  structure,  namely  that  of  the 
homogeneous  Ricatti  equation  [16].  For  a  very  long  time  (at  least 
over  100  years)  it  has  been  known  that  equations  of  this  type  may 
be  transformed  to  a  first  order  linear  form  by  introducing  a  new 
unknown  equal  to  the  reciprocal  of  the  original.  The  general 
solution  can  therefore  be  easily  written  down  as 

r 


°'1  *  M-f  i  5F  dr')fC  •  2/l ?(-/  E  3F  ir")dr' 


(8) 


where  C  is  the  integration  constant.  If  k  is  regarded  as  constant, 
which  appears  not  to  be  unduly  restrictive  and  shall  be  assumed 
here  for  brevity,  a  well  known  form  of  the  solution  may  be  obtained: 

exp[(S  -  S  )/k] 

*  *  7 - ? - 2 - y 

K1 "  i/ex pt<s  -  so)/k}dr,j 

ro 

where  <Jg  *  c(rQ) .  The  first  appearance  of  Eq.  (9)  or  its  equivalent 
in  the  literature  on  remote  sensing  was  apparently  in  1954  in  the 
context  of  rain  intensity  measuresieats  by  radar  at  attenuating 
wavelengths  [17]. 
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It  has  since  re-emerged  in  several  articles  on  the  interpretation 
of  lidar  measurements  [3,  18-20]. 

In  spite  of  the  evident  theoretical  superiority  of  Eq.  (9) 
over  the  slope  method  (which  corresponds  to  setting  k  =  0  in 
Eq.  (7)),  it  is  the  latter  method  which  is  most  often  used.  This  is 
because  Eq.  (9)  has  a  tendency  to  produce  at  best  marginal  results, 
and  in  practice  has  likely  been  more  a  source  of  frustration  than  a 
useful  tool  for  analyzing  radar  or  lidar  returns.  For  example,  in 
their  1954  article  referred  to  above,  Hitschfeld  and  Bordan  [17] 
concluded  it  was  probably  not  possible  to  calibrate  a  radar  set 
accurately  enough  to  make  use  of  the  solution,  and  that  rainfall 
measurements  made  without  correcting  for  attenuation  via  the  solu¬ 
tion  are  in  many  cases  more  accurate  than  the  corrected  values. 

Worse  yet,  others  have  noted  the  solution  may  lead  to  ”...  absurdly 
large,  infinite,  or  negative  values...”  [18]  and  "...  physically 
meaningless...”  [21]  results.  Others  have  avoided  such  behavior 
only  by  using  unrealistically  large  values  of  k  [3]. 

There  is  surprisingly  little  comment  in  the  literature  on  the 
reasons  for  the  failure  of  Eq.  (9).  It  seems  only  to  be  somewhat 
vaguely  attributed  to  the  omission  of  multiple  scattering  effects 
[1,  7].  However,  since  the  slope  method  suffers  from  the  same 
deficiency,  but  with  apparently  much  less  drastic  consequences, 
this  explanation  is  not  very  convincing. 


-rr 


.;,3r 
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Unfortunately,  only  a  few  relevant  studies  on  the  possible 
importance  of  multiple  scattering  are  available.  In  one  of  these, 
Viezee  et.  al.  [7]  compared  lidar  and  transmissometer  measurements 
in  dense  fog  and  found  an  apparent  10  to  45%  overprediction  of 
lidar-derived  visibilities  using  the  slope  method.  They  conjec¬ 
tured  this  discrepancy  was  due  to  the  influence  of  forward  and 
multiple  scattering  and  proposed  an  empirical  correction  to  the 
slope  method  for  use  under  turbid  atmosphere  conditions.  On  the 
other  hand,  they  also  noted  that  available  theoretical  descriptions 
of  multiple  scattering  [22-24]  could  not  account  for  the  observed 
discrepancies  in  the  lidar  and  transmissometer  data.  A  later  Monte 
Carlo  simulation  of  second  and  third  order  multiple  scattering  in 
dense,  homogeneous  fog  led  to  the  conclusion  that  multiply  scattered 
radiation  will  cause  the  slope  method  to  be  in  error  by  less  than 
about  10%  for  visibilites  of  the  order  of  100m  [25]. 

From  these  studies  it  appears  unlikely  that  even  for  a  dense 
dispersion  the  contribution  of  multiply  scattered  radiation  could 
make  a  crucial  difference  in  the  applicability  of  Eqs.  (1)  or  (9). 
Therefore,  although  it  would  certainly  be  desirable  to  replace 
Eq.  (1)  with  a  new  governing  form  containing  higher  scattering 
approximations  (perhaps,  for  example,  along  the  lines  recently 
outlined  by  Samokhvalov  [26]),  there  seems  at  present  no  justifi¬ 
cation  for  regarding  the  inclusion  of  multiply  scattered  radiation 
effects  as  the  sine  qua  non  for  the  inversion  of  lidar  signals  from 
a  markedly  inhomogenous  atmosphere.  (In  this  regard  it  should  be 
recalled  also  that  the  traditional  description  of  Eq.  (1)  as  the 
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single  scattering  lidar  equation  is  a  misnomer,  since  Eq.  (1) 
incorporates  the  assumption  that  backscattered  photons  are  also 
attenuated;  hence  it  does  include  some  multiple  scattering  effects.) 

From  a  purely  mathematical  point  of  view,  it  is  easy  to  see 
the  problem  with  Eq.  (9):  Since  on  average  the  signal  decays  with 
range  beyond  rQ  due  to  attenuation,  o  is  determined  as  the  ratio  of 
two  numbers  which  each  become  progressively  smaller  with  increasing 
r;  furthermore,  the  denominator,  which  must  approach  zero  at  nearly 
the  same  rate  as  the  numerator,  is  expressed  as  the  difference 
between  two  relatively  large  numbers.  Such  structure  produces  a 
strong  tendency  for  instability  and  suggests  that  unattainable 
accuracy  in  the  determination  of  Oq  may  often  be  requisite  for 
avoiding  a  singularity,  even  for  signals  which  are  free  of  noise. 

The  above  description  may  be  illustrated  quantitatively  by 
considering  the  growth  of  a  small  perturbation  in  0  due  to  an  error 
do  in.  the  determination  of  ff^.  For  the  same  signal  let  0  be  the 
solution  corresponding  to  cq,  and  a'  ~  a  +  5  be  the  solution  cor¬ 
responding  to  Oq  =  0Q  +  6q.  Then  from  the  integrated  form  of 
Eq.  (7)  it  follows  that 


which  implies 


(1  +  5/a)  =  (1  +  <5o/ao)  exp  (| 


r 

o 


(ID 


For  simplicity  consider  a  homogeneous  atmosphere  with  a  =  a  .  Then 
by  differentiating  Eq.  (11)  one  obtains 


&  =  2 
dr  k 


%C(C  - 


1) 


where  £  =  (1  +  S/a^).  This  is  a  homogeneous  Ricatti  equation,  like 
Eq.  (7),  with  a  solution  given  by 


(1  -  I-)’1  = 

o 


exp 


2a (r  -  r  ) 
0  o 

k 


(12) 


From  this  expression  it  can  be  seen  that  an  underestimate  of  a 
(5q<0)  leads  to  5-*~cro  (i.e.,  <T*i I)  as  r-*».  On  the  other  hand,  if  CTq 
is  overestimated  (5o>0),  then  <?■*»  within  a  finite  distance  given  by 


hr 


(13) 


-1  “2 

For  example,  if  aQ  =  10  km  ,  k  =  1,  and  50/^0  *  10  »  then  the 
solution  has  a  singularity  within  about  the  next  231  meters;  also, 
for  r>rQ  +  hr  the  solution  is  negative.  This  example  is  also  shown 
in  Figure  1,  where  the  unit  of  length  for  range  plotted  along  the 
abscissa  has  been  set  equal  to  0.010  km  (10m),  and  attenuation  per 
km  is  plotted  along  the  ordinate.  (This  same  choice  .  ?  units  is 
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used  for  all  the  theoretical  curves  of  c  =  a(r)  given  in  this 

paper.  Corresponding  curves  of  S  -  Sq  vs.  r  use  the  same  range 

scale;  S  -  Sq  is  of  course  dimensionless.)  The  tendencies  shown  in 

Figure  1  are  accentuated  by  larger  aQ/k  (lower  visibilities)  and 

larger  6  /a  (poorer  estimates  of  a  ). 
o  o  o 

Finally,  substitution  of  Eq.  (12)  back  into  Eq.  (10)  reproduces 
the  original  signal,  S  -  Sq  =  -2aQ(r-ro),  independently  of  the 
value  of  6q.  Therefore,  two  main  points  should  be  emphasized 
regarding  these  results:  1)  Eq.  (9)  is  -  to  loosely  paraphrase  a 
terminology  used  in  analogous,  though  generally  more  complicated 
circumstances  -  "ill  constructed",  in  that  small  differences  in  the 
choice  of  boundary  value  <J0  provide  no  assurance  that  the  corre¬ 
sponding  solutions  will  remain  close  for  r>rQ.  2)  Closeness  of 
the  S(a)  curve,  reconstructed  from  the  solution  for  a,  to  the 
original  S  curve  is  insufficient  to  guarantee  the  reasonableness  of 
the  solution.  (Such  closeness  has  been  used  in  the  past  as  a  test 
of  validity  of  the  solution  {3].)  Because  of  this  behavior  one 
would  expect,  and  experience  has  shown, -that  Eq.  (9)  by  itself  is 
of  very  little  practical  value. 

III.  A  New  Solution  Form 

It  is  fortunately  quite  easy  to  select  a  different  and  more 

appropriate  solution  form  than  Eq.  (9).  One  merely  has  to  evaluate 

the  integration  constant  C  in  Eq.  (8)  in  terms  of  a  reference  range 

r  such  that  the  solution  is  generated  for  rfr  ,  rather  than  for 

rir  as  before.  For  constant  k  the  result  is 
o 
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a(r) 


exp[(S  -  Sm)/k] 


(14) 


where  S  =  S(r  )  and  a  =  a(r  ).  This  seemingly  innocuous  change 
m  m  m  m 

from  Eq.  (9)  makes  a  very  significant  difference  in  the  behavior  of 

the  solution.  As  r  decreases  from  r  ,  a  is  now  determined  as  the 

m 

ratio  of  two  numbers  which  each  become  progressively  larger,  so 
that  stability  and  accuracy  are  easy  to  maintain.  The  form  of  the 
denominator  also  indicates  that  the  dependence  of  the  solution  on 
ct^  decreases  with  decreasing  r. 


The  contrasting  behavior  of  Eqs.  (9)  and  (14)  is  illustrated  in 
Figure  2.  (In  this  and  several  subsequent  figures,  displays  are 
given  of  various  "inversions"  (solutions  of  Eqs.  (9)  and  (14))  of 
signals  generated  by  Eq.  (10)  in  response  to  specified  a  distribu¬ 
tions.  The  value  k  =  1  was  used  in  the  computations,  except  where 
otherwise  indicated.  (The  choice  of  k  is  of  course  not  important 
so  long  as  the  same  value  is  used  for  generating  the  signal  as  for 
inverting  it.))  Figure  2(a)  shows  the  signal  response  to  the  plat¬ 
form  distribution  of  a  given  in  Figure  2(b).  Also  in  Figure  2(b) 
are  shown  the  signal  inversions  from  Eq.  (9)  for  boundary  values 
which  are  in  error  by  ±1%.  Figure  2(c)  displays  the  corresponding 
inversions  from  Eq.  (14)  for  boundary  values  which  are  in  error  by 
±50%.  The  relatively  small  effect  of  a  poor  boundary  value  estimate 
on  Eq.  (14)  is  obvious. 
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Analogous  differences  in  the  capacities  of  the  inversions  to 
survive  simulated  signal  noise  are  shown  in  Figure  3.  Figure  3(a) 
displays  three  superposed  signals,  the  first  being  due  to  a  constant 
attenuation  of  10  km  *,  and  the  others  being  like  the  first  except 
for  ±10%  "blips"  on  the  interval  (30,32).  Figures  3(b)  and  3(c) 
give  the  respective  corresponding  signal  inversions  from  Eqs.  (9) 
and  (14),  wherein  the  correct  boundary  value  has  been  used  for  all 
solutions.  It  can  be  seen  that  the  solution  of  Eq.  (9)  is  obliter¬ 
ated  for  ranges  beyond  the  point  noise  is  encountered,  whereas 
Eq.  (14)  displays  a  strong  tendency  to  recover  from  signal  errors. 

The  effect  of  an  incorrect  value  of  k  is  illustrated  in  Figure 
4,  where  again  the  computations  are  based  on  the  platform-shaped 
distribution  of  attenuation  shown  in  the  prior  figures,  and  k  s  1 
was  used  to  generate  the  signal  (shown  in  Figure  2(a)).  In  Figure 
4(a)  it  can  be  seen  that  as  soon  as  the  signal  slope  varies,  so 
that  the  value  of  k  enters  into  the  calculations,  the  inversion 
based  on  Eq.  (9)  fails.  The  much  weaker  impact  on  Eq.  (14),  illus¬ 
trated  in  Figure  4(b),  indicates  great  accuracy  in  the  determination 
of  k  is  not  required. 

In  summary,  these  examples  show  that  Eq.  (14)  is  relatively 
insensitive  to  the  kinds  of  errors  that  are  likely  to  effect  the 
inversion  of  real  signals.  Especially  encouraging  is  the  tendency 
of  Eq.  (14)  to  approach  the  correct  solution  curve  in  spite  of  a 
poor  estimate  of  the  boundary  value,  a  .  This  raises  the  hope  that 
difficult,  accurate  "lidar  calibrations"  or  independent  measurements 
of  some  optical  parameter  at  a  reference  point,  or  through  a  given 
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layer,  may  not  be  necessary,  at  least  for  the  majority  of  applica¬ 
tions.  In  the  next  section,  the  question  of  how  to  make  a  reasonable, 
self-contained  estimate  of  0^  from  the  signal  alone  is  considered. 

IV.  Estimation  of  a  :  Generalization  of  the  Slope  Method 

m 


It  would  appear  to  be  a  relatively  straightforward  matter  to 

obtain  a  good  estimate  for  a  ,  in  view  of  the  following  circumstance: 

m 

Assuming  the  validity  of  the  lidar  equation,  Eq.  (1),  and  the 

constitutive  relation,  Eq.  (6),  and  assuming  also  that  a  varies 

linearly  over  a  specified  interval  (r  ,r.  ),  then  it  is  possible  to 

express  a  solely  in  terms  of  the  signal  over  the  interval.  This 

follows  directly  from  integration  of  Eqs.  (9)  and  (14)  over  (r  ,r.), 

with  S  -»S  =  S(r  )  and  S  -»S,  =  S(r.  ): 
o  a  a  mb  b 


r  b  -k  f  ^°a 

J  adr  =  £n  1  -  J  exp[ (S  -  S^/kJdr'  ,  (15) 


j '  'crdr  =  |  in  £l  +  J  exp  [(S  -  S^/kJdr'J  , 


(16) 


where  a  =  a(r  )  and  a,  -  a(r.  ).  Because  of  the  assumed  linear 
a  a  b  b' 

variation  of  a  over  (  rg ,  r^ ) ,  the  average  value  of  a  on  the  interval 
is  just  a  =  (oa  +  /2 .  Therefore,  Eqs.  (15)  and  (16)  may  be 

combined  to  predict  the  values  of  a,  or  For  exasq>le,  a  is 
obtained  from  the  solution  of  the  equation 


q  =  (1  -  exp(-fl))  (exp(fl)  -  1) 
2Iab  2Iba  ’ 


(17) 
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where 


ft  = 


2a<rb  _  ra} 


(18) 


Xab  =  (rb  '  ra)_1  f  exTP[(S  '  SJ/kldr’, 


(19) 


and 


- 1  r  D 

Xba  =  ^rb  '  ra^  J  exP^S  "  sb)/kldr*  =  Iafa  exp[(Sa  -  Sb)/k]  (20) 


Since  the  assumption  that  a  is  linear  will  become  better  with 
decreasing  interval  size,  the  application  of  Eqs.  (17)  -  (20)  over 
a  succession  of  small  intervals  would  appear  in  principle  to  consti¬ 
tute  an  inversion  of  the  lidar  signal  which  does  not  require  any 
information  beyond  that  contained  in  the  signal  itself.  In  practice 
however,  the  local  structure  of  the  signal  is  not  known  well  enough 
to  ensure  the  success  of  such  a  method.  This  can  be  seen  by  consid¬ 
ering  the  form  of  the  solution  to  Eqs.  (17)  -  (20)  for  the  case 
that  ft«l,  i.e.,  for  intervals  Ar<<k/2o.  By  expansion  of  1^  and 
I^a  to  include  terms  proportional  to  Ar2,  the  solution  for  ff  is 
found  to  be 


.  -3(S'  +  S'.  ) 

a  =  - * - L. 
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1  + 


1 - 1* -  ( s'2  +  s;2  +  k  (s"  +  sn) 

L  9(Sa  +  S£)2  '  *  “  *  **  1 J 


,(21) 


where  S'  =  (dS/dr)r  ,  SJ  =  (dS/dr)r.  ,  S"  *  (d2S/dr2)r  ,  and 

S  D  I  I 

&£  =  (d2S/dr2)r^.  This  generalization  of  the  slope  method  result, 
Eq.  (5),  is  certainly  more  rigorous  in  its  account  of  the  local 
geometry  of  the  signal.  Unfortunately  however,  the  new  terms 
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T 


V 


ffikubi'.  * 


representing  signal  curvature  are  extremely  difficult  to  estimate, 
so  that  point-by-point  application  of  Eqs.  (17)  -  (20)  or  (21)  can 
generally  be  expected  to  provide  little  real  improvement  over  the 
slope  method. 

An  example  of  the  pomc-by -point  use  of  Eq.  (21)  is  shown  in 
Figure  5(b).  For  tho  platform  distribution  used  to  generate  the 
lidar  signal,  the  curvature  is  zero  everywhere  except  at  the  points 
where  the  slope  changes  abruptly,  and  at  these  points  it  becomes 
infinite.  The  obvious  departure >  of  the  inversion  curve  shape  from 
the  platform  distribution  are  due  to  the  effective  numerical  diffusion 
of  the  input  delta  function  distribution  of  curvature.  Thus  even 
for  this  theoretical  example  wherein  the  lidar  signal  appears  quite 
smooth  and  is  known  to  a  high  degree  of  accuracy,  significant 
errors  occur  in  the  estimation  of  signal  curvature. 

However,  in  defense  of  Eq.  (21)  it  is  also  worth  noting  that 
its  point-by-point  application,  in  conjunction  with  an  algorithm  to 
smooth  the  signal,  will  generally  permit  a  good  recovery  of  the 
broad  features  of  an  attenuation  distribution.  An  example  to  this 
effect  is  given  in  Figure  5(c).  For  the  inversion  shown  the  signal 
in  Figure  5(a)  was  smoothed  15  times  (according  to  the  simple  scheme 
that  Si-*(Si+1+Si+S^_j)/3) ,  subject  of  course  to  the  constraint  of 
fixed  curve  end  points.  In  general  it  appears  adequate  to  smooth 


the  signal  until  |s"/S'2|=10  1.  Finally,  the  inversion  displayed  in 
Figure  5(c)  also  produces  nearly  the  same  a  for  the  entire  range  as 
does  the  input  distribution.  This  happens  because  a  depends  only  on 


situations  in  which  the  endpoint  values  of  attenuation  are  known  to 
be  approximately  equal  (e.g.,  a  localized  obscurant  in  otherwise 
clear  air),  the  use  of  Eq.  (21)  with  signal  smoothing  will  not 
result  in  a  significant  error  in  the  estimate  of  total  optical  depth, 

rb 

J  adr  =  a(rb-ra) . 


Although  Eqs.  (17)  -  (20)  or  Eq.  (21)  can  be  used  to  provide 
estimates  of  a^,  the  above  discussion  suggests  there  may  on  occasion 
be  instability  problems  for  signals  with  "rough"  structure.  Greater 
stability  may  of  course  be  achieved  by  choosing  a  larger  integration 
interval,  but  offsetting  this  is  the  decreasing  likelihood  that  a 
remains  linear  over  larger  intervals.  An  alternative  which  is  not 
as  impaired  by  the  attempt  to  incorporate  more  signal  information 
over  a  larger  interval  is  the  following:  From  Eq.  (16)  one  can 
write 

r  r  r  Tk_e*p(Vk)+  /Cexp(S/k)dr'  ' 

f  b<Jdr  =  f  adr  -  f  ^dr  =  \  £n  - — - — -  (22) 

j  rJ  /  k  e*p(S  /k)  +  rc 

ra  ra  rb  - r— S -  +  /  exp(S/k)dr' 


where  rc>rij>ra'  *or  rc>>ra’  rb  tbe  hand  side  of  this  expres¬ 

sion  will  be  only  weakly  dependent  on  a^.  Imagine  now  a  hypothet¬ 
ical  linear  extension  of  the  signal  curve  beyond  rB  with  a  slope  -A 
equal,  for  example,  to  the  mean  slope  of  the  curve  over  the  range 

of  interest  for  r  SrSr  ;  i.e.,  let  r  >r  >r.  and  set  S  =  8  -  A(r  -  r  ) 

on  c  a  d  a  a 

for  r  Sr<r  .  Then 

m  c 


so  that  for  r  >>r  Eq.  (22)  becomes 
cm 


r 

/m 

■ 


exp[(S  -  Sm)/k]dr’ 


+  f“exp[(S  -  Sffl)/k]dr' 
rb 


(23) 


By  substituting  this  result  into  Eq.  (16),  a  corresponding  estimate 
of  may  be  obtained. 


Trial  applications  of  Eq.  (23)  have  shown  that  fairly  good 
results  may  usually  be  achieved  by  setting  A  =  (SQ  -  Sffl)/(ro  -  rQ) 
and  choosing  rb  -  rQ  to  be  as  large  as  within  a  few  percent  of 
r  -  r  .  The  interval  (r  ,r,  )  may,  if  desired,  be  taken  as  small 

no  <3.  o 

as  the  basic  increment  between  data  points.  One  may  then  use  o^  so 

obtained  along  with  Eq.  (14)  (letting  r  r^,  etc.)  to  generate  the 

solution  over  (rQ,rb).  Alternatively,  one  may  further  approximate 

o=o,  and  generate  the  solution  over  (r  ,r  ). 
mb  °  o  m 


V.  Inversion  Examples  and  Discussion 


Some  examples  of  inversions  generated  according  to  the  procedure 
described  above  are  shown  in  Figures  6-8.  In  terms  of  the  range 
variable  x  *  r/10  (r  in  meters)  plotted  along  the  abscissa,  the 
inversions  shown  were  obtained  from  Eqs.  (14),  (16),  and  (23)  using 
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x  *  90,  x.  =  91,  and  x  =  100.  (However,  it  should  be  ooted  that 
a  o  m 

these  choices  are  not  special  ones  required  to  produce  the  degree 

of  accuracy  displayed  in  the  figures;  other  similar  choices  produce 

similar  results,  and  indeed  such  relative  insensitivity  to  the 

choices  of  x^  and  x^  provides  a  qualitative  indication  that  the 

inversions  are  not  seriously  in  error.)  The  overall  accuracy  of 

the  inversions  can  be  seen  to  be  quite  good.  For  the  distribution 

shown  in  Figure  6(b)  the  average  value  of  the  input  attenuation 

over  the  entire  range  is  =  32.5  km  \  whereas  the  inversion 

value  is  ^out  =  32.8  km  * .  Also  shown  in  Figure  6(b)  is  an  inversion 

based  on  the  smoothed  signal  in  Figure  6(a)  (smoothed  10  times  in 

accordance  with  the  prescription  given  earlier) ;  for  it  the  average 

attenuation  is  a  .  =  33.2  km  1.  For  the  distribution  shown  in 
our 

Figure  7(b)  the  corresponding  values  are  a ^  =  38.5  km  1  and 
aQut  =  38.8  km  \  the  difference  in  this  case  being  due  almost 
entirely  to  a  slight  inversion  misrepresentation  at  the  bottoms  of 
the  ''troughs'*. 

It  is  clear  from  Figures  6  and  7  that  the  inversion  procedure 
can  recover  considerable  detail  from  the  signal.  In  this  regard  it 
is  also  noteworthy  that  special  numerical  processing  is  not  necessary 
to  obtain  good  resolution.  The  computations  used  to  generate  the 
figures  involved  only  one  slightly  nonstandard  procedure,  namely 
that  of  fitting  the  signal  curve  point-by-point  with  a  composite 
cubic  spline  [26],  in  order  to  provide  greater  accuracy  in  the 
interpolation  of  the  signal  between  base  points.  The  subsequent 
numerical  Integrations  were  all  carried  out  via  Simpson's  rule  in 
single  precision.  However,  nearly  equivalent  accuracy  can  be 
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obtained  by  using  just  single  precision,  point-by-point  trapezoidal 
rule  integration  without  interpolation.  This  simpler  approach  should 
be  more  than  adequate  in  applications  in  view  of  the  many  other 
relevant  theoretical  and  experimental  limitations. 

Examples  of  inversions  under  high  visibility  conditions  are 

shown  in  Figure  8.  For  this  case  with  a  S  10  *  km  1 ,  the  structure 

visible  in  the  signal  is  due  to  the  predominance  of  the  first  term 

on  the  right  hand  side  of  Eq.  (7)  or  (10);  i.e.,  the  signal  variation 

is  caused  mainly  by  changes  in  the  fractional  gradient  of  attenuation, 

rather  than  by  changes  in  attenuation  magnitude.  Consequently,  the 

slope  method  of  inversion  applied  to  small  intervals  gives  erroneous 

results  for  a  situation  ironically  representative  of  those  for 

which  it  is  generally  assumed  most  applicable.  The  failure  of  the 

point-by-point  application  of  the  slope  method  is  shown  vividly  by 

the  dashed  curve  in  Figure  8(a).  Even  if  one  were  to  ignore  signal 

curve  portions  where  negative  values  of  a  are  predicted  (a  strategy 

used  in  the  past  by  slope  method  practitioners),  the  predicted 

positive  excursions  in  cr  are  grossly  overstated.  (On  the  other 

hand,  it  should  also  be  noted  that  the  slope  method  applied  to  the 

entire  range  interval  provides  an  excellent  estimate  for  a.)  The 

fine  structure  of  the  input  and  inversion  distributions  is  revealed 

in  the  expanded  ordinate  scale  of  Figure  8(b).  The  corresponding 

values  of  average  attenuation  are  O.  -  0.101  km  *  and  a  .  =  0.102  km  1 

in  out 

An  interesting  new  feature  shown  here  is  the  tendency  of  the  misfit 

at  the  boundary  point  r  to  carry  on  through  the  length  of  the 

n 

inversion.  There  is  insufficient  optical  depth  for  the  inversion 
method  to  rapidly  approach  the  true  solution  curve. 
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Figures  9  and  10  illustrate  situations  where  an  additional 
small  computational  step  can  improve  the  accuracy  of  the  inversion. 

If  o  increases  or  decreases  significantly  over  the  range  (r  ,r  ), 

o  m 

then  the  use  of  Eq.  (23)  as  described  earlier  biases  the  estimate 

of  a  too  much  toward  the  value  of  a  over  (r  ,r  ).  Thus  in  Figure 
m  o  m 

9(b)  a  is  too  small,  whereas  in  Figure  10(b)  it  is  too  large.  In 
m 

such  situations  where  the  inversion  reveals  a  fairly  systematic 
trend  in  a  (or  if  such  behavior  is  known  independently,  or  can, 
e.g.,  be  discerned  directly  by  simple  inspection  of  the  signal),  a 
better  estimate  of  cr  can  be  obtained  by  iterating  Eq.  (23)  with  A 
replaced  by  twice  the  current  iterative  value  of  ct^  (cf.  Eq.  (5)). 
This  gives  greater  weight  to  the  signal  information  available  near 
r  .  An  indication  of  the  improvement  this  strategy  can  bring  about 

TQ 

is  illustrated  in  Figures  9(b)  and  10(b).  A  simpler  alternative 

which  is  appropriate  when  it  is  known  that  o  is  in  a  separate 

regime  on  some  interval  (ra>r(n)  near  r^  (e.g.,  clear  air  on  the  far 

side  of  a  smoke  cloud,  or  clear  air  above  an  inversion)  is  simply 

to  set  A  =  (S  -  S  )/(r  -  r  )  in  Eq.  (23). 

a  m  m  a 

As  a  final  example,  inversions  of  a  real  lidar  return  from  fog 
are  shown  in  Figure  11.  The  laser  used  emitted  pulses  averaging  10 
miliijoules  in  6  nanoseconds  at  1.06  pm,  and  a  20  megahertz  sampling 
rate  transient  recorder  was  used  to  produce  sample  points  spaced 
7.5  m  apart  over  the  lidar  return  [27).  The  initial  increase  in 
the  lidar  signal  shown  in  Figure  11(a)  is  due  to  increasing  (7,  and 
is  not  caused  by  incomplete  overlap  of  the  transmitter  and  receiver 
fields  of  view.  The  inversions  are  based  on  Eqs.  (14),  (16),  and 
(23),  with  A  =  (S15  -  S30)/(r3o  “  c\^-  In  Fi«ure  H(b)  are  shown 
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the  results  for  k  =  1.0  and  0.67;  the  former  value  is  probably 
better  for  fog,  but  in  any  case  the  inversions,  as  demonstrated 
earlier,  do  not  depend  strongly  on  the  choice  of  k.  For  either 
value  of  k  the  average  attenuation  is  a  =  13.0  km  1 .  This  result 
can  be  compared  to  the  visibility  as  measured  by  a  transmis some ter 
during  the  same  experiment.  The  transmissometer  visibility  v  (km), 
based  on  a  constrast  threshold  of  0.05  so  that  v  =  3.0/a,  is  0.20  km 
[27].  The  corresponding  value  from  the  lidar  inversion  is  v  =  0.23  km. 
The  extent  of  agreement  is  as  good  as  could  be  expected,  given  just 
the  uncertainties  associated  with  the  experimental  data. 

The  support  and  encouragement  of  W.  J.  Lentz  and  J.  S.  Randhawa 
are  gratefully  acknowledged.  This  work  was  performed  under  contract 
to  the  U.  S.  Army  Atmospheric  Science  Laboratory,  White  Sands  Missile 
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3  Effect  on  inversions  of  simulated  s 
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Figure  4  Effect  on  inversions  of  value  of  k. 
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Figure  5.  (Gont) 


38 


o 

_ T _ T _ f ■ _ T _ ■ 

^0.00  2aoo  4aoo  so.oo  sooo  100.00  120.00 

RANGE  (PER  10m) 

Figure  6  Inversions  based  on  Eqs.  (14),  (16),  and  (23)  for  an  unsmootbed 
and  smoothed  signal  representing  a  strongly  inhomogeneous 
atmosphere. 
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Figure  7  Inversions  based  on  Eqs.  (14),  (16),  and  (23)  for  a  very 
dense,  strongly  inhomogeneous  atmosphere. 
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Figure  8  Inversions  based  on  Eqs.  (14),  (16),  and  (23),  and  on  the 
slope  method*  for  a  case  of  high  visibility 
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